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ABSTRACT 



Context. Based on the viscous disk theory, a number of recent studies have suggested there is large scale meridional circulation in 
protoplanetary disks. Such a flow could account for the presence of crystalline silicates, including calcium- and aluminum-rich 
inclusions (CAIs), at large distances from the sun. 

Aims. This paper aims at examining whether such large-scale flows exist in turbulent protoplanetary disks. 
Methods. High-resolution global hydrodynamical and magnetohydrodynamical (MHD) numerical simulations of turbulent pro- 
toplanetary disks were used to infer the properties of the flow in such disks. 

Results. By performing hydrodynamic simulations using explicit viscosity, we demonstrate that our numerical setup does not 
suffer from any numerical artifact. The aforementioned meridional circulation is easily recovered in viscous and laminar disks 
and is quickly established. In MHD simulations, the magnetorotational instability drives turbulence in the disks. Averaging 
out the turbulent fluctuations on a long timescale, the results fail to show any large-scale meridional circulation. A detailed 
analysis of the simulations show that this lack of meridional circulation is due to the turbulent stress tensor having a vertical 
profile different from the viscous stress tensor. A simple model is provided that successfully accounts for the structure of the 
flow in the bulk of the disk. In addition to those results, possible deviations from standard vertically averaged a disk models 
are suggested by the simulations and should be the focus of future work. 

Conclusions. Global MHD numerical simulations of fully ionized and turbulent protoplanetary disks are not consistent with the 
existence of a large-scale meridional flow. As a consequence, the presence of crystalline silicates at large distance for the central 
star cannot be accounted for by that process as suggested by recent models based on viscous disk theory. 



1. Introduction 

In the past few years, more and more evidence has accu- 
mulated that suggests there are crystalline solid particles 
in the outer regions of protoplanetary disks. Spitzcr ob- 
servations in the 5-35 /jm spectral range have revealed 
crystalline signatures i n a large fraction of T-Tauri stars 
dBouwman et al.l l2~008; Olo fsson et al. I l2009t I Sargent et al 



20091) . Comets, which are believed to form far away from 



the Sun, are also known to show high crystallin ity values 
(jCrovisier et al.lll997tlWooden et al.lll999ll2007t ). This has 
recently been independently confi rmed by the returned 
samples of the Star dust mission ( Brownlee et al. 2006; 



Zolenskv et al Mcteoritic records also indicate that 



calcium- and aluminum-rich inclusions (CAIs) are a com- 
mon component of chondrites collected from parent bod- 
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ies thought to originate in the main asteroid belt. CAIs, 
as well as the crystalline silicates observed by the Spitzer 
telescope, are believed to have formed in the inner parts of 
protoplanetary disks (R < 1 AU) where temperatures are 
in excess of ~1000 K as required for their formation from 
amorphous silicates. This general trend of finding crys- 
talline silicates at large distances from the central star re- 
quires a mechanism able to transport them from the disk's 
inner parts to their outer regions. 

Several proc esses have been s ugges ted to account for 
such transport. IShu et al.l (|1996l 120011 ) have invoked the 
action of powerful winds from the youn g stars, the so- 
called X-wind model (see also Hu| 2010h . Many authors 
have, on the other hand, relied on processes operating 
within the bulk of the disk itself. Indeed, the flow in ac- 
cretion disks is believed to be turbulent in order for an- 
gular momentum to be efficiently transported outward. 
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The turbulent nature of the flow results in an efficient 
diffusion of small solids in the disk, potentially transport- 
ing some of them to large heliocentric distances. Using 
the standard ID (i.e . verti c ally integrated) a disk mode l 
(jShakura fe Sunvaevl [l973t iLvnden-Bell fe Pringld Il974) 
as a theor e tical basis to model the effect of the t urbu- 
lence. iGaill (|200ll ) and lBockelee-Morvan et al.l (|2002h have 
quantified the resulting concentration of crystalline sil- 
icates in the outer parts of protoplanetary disks. The 
conclusion is that turbulence alone appears unable to 
transport enough solid material out to large distances 
because the inward flow associated with mass accretion 
onto the central object dominates over a long time and 
reduces its effect. As a remedy, several papers have in- 
voked the existence of a large-scale meridional flow in 
protoplanetary disks. Such a meridional circulation nat- 
urally arises out of viscous disk models that extend the 
standard a disk model to two dimensions, t he radial po- 
sition R and distan c e to the midplane Z (Urpin 19841 
Siemiginowska 19881: Kiev fe Linl 1992 ; Rozvczka et al 



1994 iRegev fe Gitelmanl 120021: ITakeuchi fe Linll2002l) . In 

these models, the radial velocity is positive in the disk's 
equatorial plane while gas accretion proceeds through 
the disk surfaces. This outwardly directed mass flux ad- 
vectively transports solids sedimcnted in the disk equa- 
torial plane out to large distances and thus circum- 
vents the limiting effect of mass accretion mentioned 
above. Several models based on this idea have been 
developed in the last few years (ITakeuchi fc Lin [2002 : 

al2007 



Keller fc GaiJ2004tlTscharnuter fc GaiJ2007tlCieslal2007 



20091: iHughes fc Armitagd 12010) . They have indeed been 
able to successfully reproduce the amount of crystalline 
solids found in outer protoplanetary disks. 

However, all of these models rely on the a prescription, 
a large-scale model for the turbulence. While the nature 
of the turbulence in protoplanetary disks is still somewhat 
debated, it is most likely magnetohydrodynamical (MHD) 
in nature and driven by the magnetorq tational instabil- 



ity (MRI, iBalbus fc Hawlevl Il99ll Il998h . Thanks to the 



large increase in computational resources, it is now possi- 
ble to perf orm global numerical simulations of protoplane- 
tary disks ( Papaloizou fc Nelsonll2003t IFromang fc Nelson 



20061: iLvra et al.ll2008l) and to apply such simulations to 



a variety of issues related to planet format i on, includ- 
ing dust dynamics (Fromang fc Nelson! 2005 : Lvra et al 



20081 Fromang fc Nelson! 20091) . plane tesimals evolution 
dNelsonl 120051: iNelson fc Gressell [2010l) . dead-zone prop- 
erties ( Dzvurkevich et al. 12010 ), and planet-disk interac- 
tion ( Nelson fc Papaloizou 12003 : Papaloizou et al. 2004 



Nelson fc Papaloizou! I200I IWinters et all 12003!) Using 
such simulations, the nature of the large-scale flow can be 
investigated from first principles without having to rely 
on ad hoc modeling of the turbulence. It is the purpose 
of the present paper to develop such dedicated numerical 
simulations in order to validate the models of meridional 
circulation presented above that are used to explain the 
presence of crystalline solids at large distances from the 
central objects in protoplanetary disks. 



The plan of the paper is as follows. In section ® we 
describe the properties of meridional circulation as it can 
be derived from 2D viscous disk theory. Section [3] presents 
a set of numerical simulations of protoplanetary disks and 
analyzes the properties of the large-scale flow. In section^ 
an additional scries of purely hydrodynamical simulations 
are used to assess the influence of the numerical setup 
and algorithm on meridional circulation. Finally, in sec- 
tion [5j the results are discussed using a simple model and 
their consequences for crystalline silicates radial mixing 
are highlighted. 

2. Large scale flow in protoplanetary disks 

In this section, we seek to derive the equations that govern 
the large-scale radial flow in protoplanetary disks. To do 
so, we consider an axisymmetric disk and derive its global 
structure. The large-scale gas flow is determined by the 
mechanism governing angular momentum transport. Two 
cases will be considered. We first examine the case of vis- 
cous disk models in section 12.31 and adopt the standard a 
prescription for the kinematic viscosity. The analysis we 
present largely follows t he eq uations derived, for exam- 
ple, bv ITakeuchi fc Linl ( 2002 ) to which the reader is re- 
ferred for further details. Since the flow in protoplanetary 
disk is most likely turbulent, we also write in section \2. 41 
the equations that govern angular momentum transport 
in such disks, highlighting the similarities and differences 
with viscous disk models. 

2.1. Definitions 

Unless otherwise stated, we consider in this paper a 
cylindrical coordinate system (R, (j>, Z) with unit vectors 
(e.R,e0,ez)- The disk model is fully specified once the 
equation of state (EOS), midplane density, and viscosity 
are chosen. For the former, we used a locally isothermal 
EOS: the sound speed is time independent and given as a 
function of the cylindrical radius R by the power law, 

R 



(1) 



where cq is the sound speed at the fidutial radius R = Rq. 
The disk midplane density p m id is similarly given by 



R 

in which po stands for the gas density at R = Rq. 



(2) 



2.2. Density and angular velocity 

The first step is to derive the spatial distribution of gas 
density and angular velocity. They are given by the equa- 
tions of force balance in the radial and vertical direction: 

^ 2 -^ GA ll.„-^=0, (3) 



(R 2 + Z 2 ) 3 / 2 p dR 
GMZ 1 dP 
(R 2 + Z 2 ) 3 / 2 ~ ~ P ~dZ 



= 0. 



(4) 
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The second equation gives the vertical profile for the den- 
sity: 



exp 



( GM 
{ c 2 



1 



VR 2 + Z 2 



where Eq. ([2]) has been used to write the disk midplane 
density. Upon expanding Eq. ((5]) to the second order in 
Z I R, one recovers the standard Gaussian vertical profile 



exp 



z- 



H 2 ) 



(0) 



where H = c s /Q,k is the disk scale height, which is the 
ratio of the sound speed to the Keplerian angular velocity 
£Ik = y/GM/R 3 . In this last definition, G and M stand 
for the gravitational constant and the central stellar mass, 
respectively . With the definition of c s given in section |2~T1 
H can be expressed as 



H 



R \ (?+3)/2 



(7) 



where H = Co/ y/GAI/Rf, is the disk scale height at R = 
Rq. Finally, Eq. (O, along with the expression for the gas 
density, gives the spatial variations of il: 



n = fl 



K 



(1 + 9) 



qR 



1 



VR 2 

2 



z 2 



+ (P + l) "ET 



1/2 



p + q + 



qZ^_ 
2H 2 



(8) 



where the last equality results from a second-order expan- 
sion in Z/R. 

These expressions for the density and the angular ve- 
locity do not depend on the viscous or turbulent nature of 
the flow. In particular, they are expected to hold in turbu- 
lent protoplanetary disks provided turbulent fluctuations 
are properly averaged out on long timescales, as well as in 
viscous accretion disks regardless of the form of the vis- 
cous stress tensor. We show in section 13.2.11 that this is 
indeed the case. 

2.3. Angular momentum conservation in viscous disks 

A widespread and convenient way of modeling the effect 
of disk turbulence is to solve the hydrodynamic equa- 
tions adding a nonzero kinematic viscosity. Since the early 
1970's, this has led to the development of a disk models in 
which the kinematic viscos i ty is assumed to be of the form 
(IShakura fc Sunvaevllil)7l iLvnden-Bell fc Pringlelll974h 

(9) 



v = ac*H . 



In this section, we examine the consequences of such a 
prescription on the disk's large-scale radial flow. Angular 
momentum conservation writes in this case as 



v «M + Vz Jz) ^ 



dz ^ R Tz * > 



(10) 



(5) T, 



where the viscous stress components Tjj" c and T|^ c are 
given by 

= pvR- 



Re, 



T. 



Zo 



pvR 



OR 
dZ 



(11) 
(12) 



In addition, we have the orderings vz ~ (H/R)vr and 
d/dZ~(H/R)d/dR. Thus, in Eq. (0, the term v z d/dZ 
is a factor (H/R) 2 smaller than the term vnd/dR and 
can be safely neglect ed in thin disks for which H<^.R 
(|Takeuchi fc Lmll2002l) . This gives 



RpvR 



dR 2 n 
dR 



d / 3 dtt 
dR \ R p »dR 



d / 3 dfl 

dZ \ R PV BZ 



(13) 



At this point, an expansion to the second order in (Z/R) 
is performed in order to proceed further in the analysis. 
After some manipulations and using Eq. ((5]) for the kine- 
matic viscosity, the r adial velocity can be expressed as 
(jTakeuchi fc Lmll2002l) 



Vr 
Co 



Ho 
Ro 



R_ 

Ro 



9+1/2 



'dp + 2q + 6 



5^ + 9 



(14) 



In the disk midplane, the radial velocity is thus posi- 
tive whenever 3p + 2q + 6 < 0, which is readily obtained 
for standard disks parameters. For example, typical values 
generally consi dered in numerical simulations are q = — 1 
and p = - 3/2 dFromang fc Nelso"n1 l2006l 120091 iLvra et al 



a p = 

2009^ iDzvurkevich et all l2010t ). and they correspond to 



outward midplane velocities. For these typical parameters, 
the second term appearing in front of the Z dependence 
is positive. This indicates that the flow direction reverses 
in the disk's upper layers to become negative. Overall, 
this flow pattern produces a meridional circulation in the 
disk in which gas flows outward in the disk midplane and 
accretes through the disk's surface layers. 

Equations ([6]), ([8]), and (fT4|) completely determine the 
flow structure in a viscous protoplanetary disk. It is es- 
tablished on a f e w dyn amical timescales, as suggested by 
Takeuchi fc Lin (120021) and demonstrated using 2D nu- 
merical simulations bv iRozvczka et ah ( 19941 ). This does 
not mean, however, that such a disk is in steady state. 
Indeed, that is only the case for particular combinations 
of the exponents p a nd q for which the ma ss flux is in- 
dependent of rad ius ( Takeuchi fc Lin 20021 ). In a more 
genera l situation, Rozvczka et al.l ( 19941 ) and Kiev fc Linl 
(|1992l ) demonstrated using 2D numerical simulation that 
a meridional circulation is quickly established, on the one 
hand, as a dynamical response to the viscous stress. On the 
other hand, the disk structure evolve viscously on much 
longer timescales. We confirm these results using our own 
simulations of viscous disks in Sect. |U 

It is important at this stage to stress that the am- 
plitude of the radial velocity predicted by Eq. ([T4l is ex- 
tremely low. Indeed, as indicated by Eq. (|T4l . the Mach 
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number of the radial velocity is about aHo/Ro. With stan- 
dard values such as a ~ 1CP 2 and Ho/Rq ~ 0.1, this 
gives vr/co ~ 10~ 3 . In contrast, years of numerical simu- 
lations of MRI-induccd MHD turbulence have taught us 
that the amplitude of the turbulent velocity fluctuations 
are around 5 to 10% of the sound speed, i.e. almost two 
orders of magnitude larger. Detecting the signature of a 
meridional circulation in turbulent disk simulations will 
thus require extremely long simulations to properly aver- 
age out the turbulent velocity fluctuations. 

2.4. Angular momentum conservation in turbulent 
disks 

In a turbulent disk, angular momentum conservation can 
be written in a form similar to Eq. (|13[) . 



RpvR 



dR 2 fl d 



dR 



dR ^ Tr * } 



r) 

dz ^ Tz * > 



(15) 



but where T R u J b and T^"' b now stand for the turbulent 
stress tensors and are respectively given in terms of the 
velocity and magnetic field turbulent fluctuation s by the 
following expressions (|Balbus fe PaDaloizoulH999h : 



T R< p = < B R B$ - pvrSvq > 
Tzcj, = < B z Bi, - pv z 5v$ > , 



(16) 
(17) 



where 5v$ stands for the fluctuations of the azimuthal 
component of the velocity. While analytical calculations 
have shown in section [2] that the spatial variations in the 
viscous stresses T^ c and Tgg c lead to large-scale merid- 
ional circulation in protoplanetary disks, no such deriva- 
tion can be performed for turbulent disks since their radial 
and vertical profiles are unknown and highly fluctuating in 
time. One therefore has to rely on numerical simulations 
instead. In the following sections, we present the results of 
such simulations and analyze the effects of the turbulent 
stresses on the large-scale flow in turbulent protoplane- 
tary disks. 

3. Numerical simulations 

3.1. Setup 

Three models of turbulent protoplanetary disks are pre- 
se nted in this paper. The s e tup is very similar to those 
of iFromang fc Nelson! (|2006l . l2009h . The MHD equations 
are s olved with the code GLOBAL (jHawlev fc Stone 
1995h using spherical coordinates (r, 9, (f>) with unit vec- 



tors (e r ,eg,e^). The resolution is set to (N r , Ng, N^) = 
(512, 256, 256). At time t = 0, the gas density and the an- 
gular velocity are initialized such that the disk is in hydro- 
static equilibrium. This is done using Eq. ([5]) for the gas 
density and Eq. © for the angular velocity. The radial and 
meridional velocities are both set to zero. Units are such 
that GM = 1, Rq = 1, po = 1, and c = 0.1. For all mod- 
els, we used q = — 1. This value gives a steeper temper- 
ature radial profile than suggested by most realistic disk 



models. It was largely chosen for computational reasons 
since it gives a constant opening angle for the disk. Taken 
together, these relations mean that H/R = Hq/Rq = 0.1 
everywhere in the disk. The computational grid covers the 
radial extent r € [1, 10], the azimuthal range <j> G [0, w/2], 
and five scale heights on both sides of the equatorial plane: 
6 G [it/2 - 0.5,7r/2 + 0.5]. The grid resolution thus cor- 
responds to more than 25 cells per scale height in the 
meridional direction. Finally, three values were considered 
for the last parameter: p = —3/2, —2, and —5/2. They re- 
spectively give a physically reasonable power law radial 
profile for the surface density with exponents —0.5, — 1, 
and —1.5. In the rest of this paper, time is measured in 
orbital periods at the grid inner edge. Using this set of 
parameters, each si mulation can be rec a st int o physical 
units as explained in Fromang fc Nelson ( 20091 ). 

To trigger the MRI at the beginning of each run, the 
disk is initially threaded by a weak purely toroidal mag- 
netic field whose strength is such that the ratio (3 between 
the thermal and magnetic pressure equals 25 everywhere. 
Random velocity fluctuations are then added to each com- 
ponent of the velocity, with an amplitude equal to 1% of 
the local sound speed. The boundary cond itions are the 
same as used by IFromang fc Nelson! (l2006h and are peri- 
odic in azimuth and outflowing in the meridional direc- 
tion. The inner radial boundary uses the "viscous out- 
flow conditions" . in which the radial velocity is set to 
a constant value consistent with a disk theory, while all 
other variables are outflowing. The outer radial boundary 
is handled through a nonturbulent buffer zone (covering 
the range R £ [9, 10]) where motions and magnetic fields 
are damped. 

With this setup, the disk quickly becomes fully turbu- 
lent because of the MRI. The idea is then to average out 
the turbulent velocity fluctuations in order to extract the 
mean flow structure of the disk. This averaging is done 
both in space (using an azimuthal average) and in time. 
However, as detailed in section 12.31 the expected ampli- 
tude of the meridional circulation is much smaller than 
the amplitude of the turbulent velocity fluctuations. As 
a result, the raw data need to be averaged over a long 
time interval to average these fluctuations out. We typi- 
cally found that about 400 dumps evenly spaced over 200 
orbits at R = 1 were needed to extract a meaningful signal 
from the data. In the middle of the computational domain, 
say at R = 5, this corresponds to about 20 orbits. This is 
much shorter than the viscous timescale (typically a few 
hundred local orbital time) over which the disk evolves 
away from the initial state. Thus we expect that the disk 
structure in the middle of the computational domain will 
remain close to its initial state described by the power-law 
functions above. This will facilitate the comparison with 
the models presented above. 

An additional complication comes, however, from the 
violence of the MRI linear phase. Indeed, the transients 
associated with it significantly affect the disk structure 
and drive it away from its initial state by the time the 
entire disk has become fully turbulent. This prevents a 
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bar symbols denote density-weighted azimuthal and ver- 
tical averages. For example, 



Fig. 2. Disk surface density for the case p = — 2 at t — 400 
(solid line) and t = 600 (dashed line). 



detailed comparison with viscous disk models. To solve 
that problem, the simulations were stopped after 300 or- 
bits. At that time, the density and velocity were reset to 
their initial values, keeping the magnetic field to its cur- 
rent value. The models were then restarted for another 
300 orbits. Time averaging was usually started about 100 
orbits after this restarting procedure to avoid any spuri- 
ous transient that might be associated with it. As shown 
in the next section, this procedure greatly helps to main- 
tain the disk close to its initial state and allows a clean 
comparison between viscous and turbulent disk models. 

3.2. Fiducial run: p=-2 

We first present a detailed analysis of the model for which 
p = — 2, before using the other two cases to demonstrate 
the robustness of the result. All the data presented in this 
section were obtained after time averaging the raw simu- 
lation results between t = 400 and t = 600 using roughly 
400 snapshots evenly spaced in time by half an orbit. 



Sv R 5v,f, 



J J p5vii5v^d(j)dZ 
J J pd(f>dZ 



1 



2ttI] 



pdvuSv^d^dZ , 
(19) 

where the last relation stands for a definition of the disk 
surface density S. 

The radial profile of a as defined above and time aver- 
aged between t = 400 and t = 600 orbits is shown in fig- 
ure [1] (left panel). In agreement with published models of 
strat i fied protoplanetary disks ( Fromang fc NelsorJ 12006 . 
20091: iDzvurkevich et alJl20ld ISorathia et alJl2010h . a is 
on the order of 0.005. It shows fairly smooth variations 
with R that indicate a systematic decrease with radius. 
The reason for such a radial decline is still unclear and 
beyond the scope of this paper (see however a detailed 
discussion of that point in section [53)1 . 

The right hand panel of figure [TJ displays the vertical 
profile of fi at radius R = 4.5, averaged azimuthally over 
the computational box and in time between t = 400 and 
t = 600. For the purpose of comparison, the figure also dis- 
plays the prediction of Eq. ((5J). There is a small systematic 
difference between the two curves because the radial sur- 
face density profile is slightly flattened at R — 4.5 over 
the course of the simulation (see fig. 0), thereby reducing 
pressure support and increasing Q. But aside from this 
small difference, there is good agreement between the two 
curves. Finally, the quasi steady state structure of the disk 
is illustrated by figure [2 which shows the radial profile of 
the surface density at t = 400 (solid line) and t = 600 
(dashed line). The two curves differ at most by 20 to 30%, 
and the difference is even much smaller in most parts of 
the disk. These small variations in the disk surface den- 
sity over the duration of the simulations confirm that the 
disk structure remains close to its initial state, as expected 
since the viscous timescale in the middle of the computa- 
tional domain is much longer than the duration of the 
simulation. It can thus be accurately described using the 
power laws we used in section [21 which makes the present 
simulations an excellent laboratory for studying whether 
meridional circulation exists in turbulent protoplanetary 
disks. This is the purpose of the following section. 



3.2.1. Flow properties 

Before considering the possibility that a meridional cir- 
culation develop in the disk, we first describe some basic 
features of the flow that will be needed for the rest of 
the analysis. The most important of them is provided by 
the parameter a, which measures the rate of angula r mo- 
mentum transport. As in iFromang fc Nelson! (120061 ). it is 
measured as a function of radius according to 



a(R) = ttRey + CKMax = 



SvrSv^ - -fcj 



(18) 



where aricy and ctMax correspond to the Reynolds and 
Maxwell stress contributions to a, respectively. The over- 



3.2.2. Meridional circulation 

The properties of the large-scale radial flow that develops 
during the simulation arc illustrated in figure[3J For the re- 
sults to be easier to interpret, only the region \Z\ < 2.5H 
has been considered in the analysis. This corresponds to 
the midplane region of the disk where most of the mass 
concentrates. According to Eq. (fl4")) . this is also the re- 
gion where meridional circulation is expected to develop. 
Velocity fluctuations above that height are so large that 
no converged mean flow could be extracted from the sim- 
ulations, even after averaging for 200 orbits. In any case, 
because of the exponential drop in density, little mass flux 
is expected in that region. The sign of the radial velocity 
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Fig. 1. Left panel: radial profile of «R cy (dotted line), ctMax {dashed line), and a (solid line), as defined according to 
Eq. (|f 8[) . for the case p = —2. Right panel: vertical profile of 51 at R = 4.5 for the same model (solid line) compared 
with the prediction of Eq. §E§ shown using the dashed line. 




Fig. 3. Left panel: Time and azimuthally averaged radial velocity for the model p = —2. Positive velocities are marked 
with white colors, while black regions correspond to negative vr. The raw simulation data have been averaged in time 
between t = 400 and t = 600. Right panel: The solid line shows the vertical profile of the radial velocity averaged in 
time between t = 400 and t = 600, in the azimuthal direction and in the radial direction between R — 3 and R = 6. 
The dashed lines show the theoretical prediction of Eq. (jT4")) for a = 1Q~ 2 and 5 x 10 -3 , respectively . The dotted line 
simply marks the zero point as a reference. 



in the (R, Z) plane is displayed in the left hand panel 
of figure [3] Despite the very long averaging period, the 
patchy structure apparent in that figure indicates that the 
data remain very noisy even after the long average we per- 
formed. As discussed above, this is due to the large-scale 
flow having an amplitude much smaller than the turbu- 
lent velocity fluctuations. Nevertheless, it is also evident 
from the figure that the numerical simulations do not show 
any obvious signature of a meridional circulation. Such a 
feature would have been characterized by a large white 
region around the disk midplanc sandwiched between two 
black regions. In fact, it appears that many black blobs are 
clustered around the disk midplane, while the disk surface 



(i.e. the region located at larger \Z\) appears whiter. This 
is confirmed by the right hand panel of that figure, in 
which an additional radial averaging has been performed 
between R = 3 and R = 6. The resulting vertical pro- 
file of vr is compared with the theoretical expectation of 
Eq. (TT4")) for two values of a, namely 5 x 10~ 3 and 10~ 2 . 
In the simulations, vr is found to be low in the vicinity 
of the disk midplane and rises in the disk surface layers. 
It is surprising to see positive velocities at all Z since it 
indicates a mean bulk motion of the disk. This will be 
largely addressed in section I5T21 where the connections be- 
tween the present simulations and vertically averaged a— 
disks arc discussed. Here, we only comment that the large 
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Fig. 4. Vertical profile of the radial velocity for the case p = —2. In both panels, the solid line is identical to that 
plotted on the right panel of figure [3] (time average over [400,600] and radial average over [3,6]). The left panel 
investigates the sensitivity of the result on the time average while keeping the same radial range: the time average 
is taken over [400,500] (dotted line) and [500,600] (dashed line). The right panel shows data obtained with the same 
time average, [400,600], but different radial range, namely R £ [2,4.5] (dotted line) and R £ [4.5,7] (dashed line). 



difference between the simulations and the predictions of 
viscous disk theory further supports the conclusion sug- 
gested above. There is no meridional circulation in turbu- 
lent protoplanetary disks in which turbulence is driven by 
the MRI. 



As mentioned above, the amplitude of the large-scale 
flow is nevertheless much smaller than the typical turbu- 
lent velocity fluctuations. This raises the question of the 
uncertainty on the radial velocity vertical profile shown in 
the right hand panel of figure [31 To investigate this point, 
we have varied the intervals over which the time and radial 
averages are performed. The result is shown in figure E] In 
both panels, the solid curve is identical to the one plotted 
in the right hand panel of figure [3] the time average is 
performed over the range [400, 600] and the radial average 
over the range [3,6]. In the left hand panel, the dotted 
and dashed lines retain the same radial range, R £ [3, 6], 
while the time interval is respectively varied over the range 
[400, 500] and [500, 600]. In the right hand panel, the time 
average is performed over the range [400,600] while the 
radial average is done over [2, 4.5] and [4.5, 7]. As a result, 
all of these additional curves arc obtained by averaging 
the simulation results cither over a smaller radial extent 
or over a shorter time interval than used to compute the 
solid curve. They are thus expected to be less converged 
than the solid line. Nevertheless, the dotted and dashed 
curves in both panels are qualitatively similar to the solid 
line (i.e. no meridional circulation is found) and show only 
moderate quantitative difference with the solid line, indi- 
cating that the results shown in figure [3] arc converged 
fairly well. 



3.2.3. Turbulent vs. viscous torques 

The results presented above demonstrate important dif- 
ferences between the large-scale flow properties of viscous 
and turbulent disks. This is due to differences between 
viscous and turbulent stresses. Indeed, as outlined in sec- 
tion 12.41 stress tensors responsible for angular momentum 
transport have no reason to be identical in viscous and in 
turbulent disks. We now compare them in detail. 

We first consider the Z(f> components of those stresses. 
Both T^ c and T^ b are respectively given by Eq. (TTJ, 
along with the a-prescription for the viscosity and 
Eq. (fl~7f . The left hand panel of figure [5] compares the 
vertical profile of T*^ rb / 'P mld and T^ c /P mtd . As for the 
right hand panel of figure^ the simulation data have been 
averaged in azimuth over the entire computational domain 
and in radius between R = 3 and R = 6. A further aver- 
age over II snapshots evenly spaced between t = 300 and 
t = 600 was also performeqj. The two stresses display a 
comparable amplitude and a similar vertical profile, which 
is ultimately related to the sign and the amplitude of the 
vertical derivative of fl. 

We now turn our attention to the R(f> components of 
the viscous and turbulent stress tensors. These are given 
by Eq. along with the a-prescription for the vis- 



Unlike the R(j> component of the turbulent stress, the 
derivation of the Zcj) component of the turbulent stress had 
not been anticipated before the simulations were performed. 
Their calculation was done using the small number of dump 
files available after their completion rather than during the 
simulations themselves as for the Tr$ components. Although 
such a procedure unfortunately results in larger fluctuations, it 
is not prohibitive because stress tensors tend to converge faster 
than other statistical diagnostics like the mean radial velocity. 
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(solid line) and T^ % f c '/ 1 P m id (dashed line), time averaged between 
600 and radially averaged between R = 3 and R = 6. The two stresses are similar both in amplitude 



Fig. 5. Left panel: Vertical profile of T^ b /P mid 
t = 400 and t 

and in their vertical profile. Right panel: Vertical profile of T^T 1 '/ 'P m id (solid line) and I P m id (dashed line), time 
averaged between t = 400 and t = 600 and radially averaged between R = 3 and R = 6. Although they are comparable 
in magnitude, the viscous and turbulent stresses display largely different vertical profiles. 
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Fig. 6. Left panel: Spatial variation of T^ rfl '/ 'P m id in the disk's meridional plane, time-averaged between t = 400 
and t = 600. Right panel: Same as the left panel, but for the quantity T^f c / P m id- A strong vertical stratification is 
apparent in the viscous stress but appears to be much weaker in the turbulent stress. 



here). First, both the viscous and turbulent R<p compo- 
nents of the stress are much larger in amplitude than the 
Zcj) components plotted in the left hand panel. Second, 
their vertical profiless are markedly different. While T^ b 
vertical variations are weak, those of T^ c are strong. 
This is because the latter traces the large vertical gra- 



cosity, and Eq. (fT6|) . respectively. Snapshots of both ten- 
sors in the disk's meridional plane, properly averaged in 
time and azimuth, are shown in figure |H1 The left hand 
panel shows T^ b / ' P m id while the right hand panel plots 
Tjil c /Pmid (as for the left panel of figure [31 only the re- 
gion within 2.5H of the equatorial plane has been plotted 
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dient of the gas density, while the former arises mainly 
because of magnetic forces. This difference is made even 
more apparent in the right hand panel of figure [5] in which 
a further radial average between R = 3 and R = 6 has 
been performed. (The range of the x-axis has been ex- 
panded up to 4.5H compared to figure [5]) In this fig- 
ure, the solid line plots the turbulent stress while the 
dashed line represents the viscous stress tensor. As sug- 
gested by figure El the viscous stress displays a vertical 
profile reminiscent of the Gaussian density profile. The 
turbulent stress, because it is dominated by the magnetic 
forces, shows a plateau for \Z\ < 2.5H before dropping 
to lower values higher in the disk. Such a vertical pro- 
file for the turbulent stress tensor has already been re- 
ported in numerous numerical simulations of MRI-driven 
MHD turbulence, using either a local (jMiller fc Stone 
2000HHirose et al.ll2006l: iFlaig et al.ll2010h or a global ap 



proachdFromang fc Nelsonl 20061 Dzvurkevich et al.ll2010t 
ISorathia et al.ll2010h . As becomes clear in section [5TTT it is 
precisely this peculiar vertical structure in the turbulent 
stress that prevents meridional circulation from develop- 



3.3. Dependence with the surface density power law 
index 

Even though the previous results suggest a significant dif- 
ference between viscous and turbulent models, the simu- 
lations presented above still display large fluctuations. To 
assess the robustness of the result, namely the absence of 
any meridional circulation in turbulent disks, we present in 
this section two additional simulations for which p = —3/2 
and p = —5/2. These two runs serve two purposes: first, 
they probe different disk radial structure. Second, and 
maybe more importantly, they test the nature of the large- 
scale flow in disks given different realizations of the tur- 
bulent flow. 

The results of these two simulations are displayed in 
figures [7] and [5] for the case p = —5/2 and p = —3/2, 
respectively. Both figures are similar to figure [3] the left 
hand panel shows the meridional distribution of vr, with 
white regions corresponding to outward motions and black 
regions to inward flow; the right hand panel plots the ver- 
tical profile of vr radially averaged between R = 3 and 
R = 6. This is compared to the predictions of viscous disk 
theory (depicted using a dashed line). Both models con- 
firm the differences between turbulent and viscous disks. 
No meridional circulation is observed, and the results are 
similar to the fiducial model p = —2. The absence of 
meridional circulation in turbulent protoplanctary disks 
thus appears quite general (i.e. independent of the disk 
structure). 

4. Numerical checks 

The MHD simulations presented in the section above sug- 
gest that meridional circulation is not present in turbulent 
protoplanctary disks. However, because of its very small 




Fig. 9. Vertical profile of the gas radial velocity ob- 
tained in two hydrodynamical viscous simulations with 
the Pencil code (red solid line) and with the JUPITER 
code (solid blue line). The two sets of curves plot the sim- 
ulations results obtained with a — 5 x 10 -3 and 10 -2 
while the dashed lines are the analytical predictions of 
Eq. (fTl| . Numerical and analytical results are in excellent 
agreement with each other. 



amplitude, one has to be careful since a number of po- 
tential numerical issues might affect this conclusion. For 
example, the grid meridional boundaries might cause wave 
reflections that could perturb the development of that 
circulation, while the limited radial extent of the com- 
putational domain perturbs its viscous evolution. Both 
problems can potentially affect the development of merid- 
ional circulation. To address these issues, we have used 
two hydrodynamic c odes, namely the Pencil CodeH and 
the code JUPITER (Ide Val-Borro et al.lf2006h . Both codes 
were used to solve the hydrodynamic equations in spheri- 
cal geometry with a prescribed kinematic viscosity. The 
implementation of the viscous force in both codes was 
checked using the test described in Appendix [21 At t = 0, 
a disk in hydrostatic equilibrium is initialized as described 
in section l2~2l using q = —I, p = —2, and H n /R = 0.1. 
The a prescription is used with the two different values 
a = 5 x I0~ 3 and I0 -2 . The size of the grid and the 
resolution in the disk's meridional plane are identical to 
the MHD runs performed with GLOBAL. However, the 
boundary conditions differ. The Pencil simulation was 
done with periodic boundaries in 9. JUPITER in turn 
used simple reflexive boundary conditions both in colat- 
itude and radius, so that the ghost zones are filled with 
copies of the neighboring zones of the active mesh, using 
a mirror symmetry. This yields a trivial Riemann problem 
at the boundary, hence a zero mass flux and a momentum 
flux that scales with the pressure at the boundary. 

The results obtained with the two codes are summa- 
rized in figure [SI in which the vertical profile of the ra- 



See http : //www.nordita. org/ software/pencil- code 




dial velocity is plotted for the two values of a quoted 
above. All curves are averaged in azimuth and in radius 
between R — 3 and R = 6. They were obtained at times 
t = 250 (Pencil Code results), t = 37.5 (Jupiter re- 
sults, case a = 5 x 10~ 3 ) and t = 57.5 (Jupiter re- 
sults, case a = 10 -2 ). For both codes and both values 
of a, the agreement between the numerical results and 
the analytical predictions is spectacular. There is only a 
tiny mismatch above 2H, which is most likely due to the 
meridional boundary conditions. This remarkable agree- 
ment shows that meridional circulation, when it exists, is 
a robust feature of the flow that is easily recovered in nu- 
merical simulations. It is neither affected by the presence 
of artificial boundaries nor strongly influenced by the de- 
tails of the numerical algorithm. In addition, these results 
are obtained at much shorter times than the disk's vis- 
cous timescale at R = 4.5. Indeed the later amounts to a 
few hundreds orbits, while the results are obtained after a 



few tens of orbital times at that radius. This confirms the 



claims of iKlev fc Lid ([1992(1 and iRozvczka et all (|1994h 



that meridional circulation is established quickly, i.e. on a 
much shorter timescale than the disk's viscous timescales. 
As a conclusion, these viscous disk numerical simulations 
give further confidence that the lack of meridional circu- 
lation reported in section [3] is a real feature of turbulent 
protoplanetary disks. 



5. Discussion 

The numerical simulations presented above have shown 
striking differences between viscous and turbulent disks. 
In this section, we first present a simple model that ac- 
counts for the numerical results and discuss the relation- 
ship between these results and ID standard a disk mod- 
els before highlighting the consequences of our findings 
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Fig. 10. The solid line shows the vertical profile of vr 
for the case p = — 2, derived as described in figure [3] It 
should be compared with the dashed line which plots the 
prediction of Eq. (j2"5)l . using a t = 5 x 10~ 3 and Hq/Rq = 
0.1. 



for large-scale radial transport of solids in protoplanetary 
disks. 



5.1. A simple model 



In section T3. 2. 31 it was shown that meridional circulation 
fails to be established in turbulent protoplanetary disks. 
As explained in section 13.2.31 this behavior, unlike that 
of viscous disks, is due to differences between the viscous 
and turbulent stress tensors that are responsible for an- 
gular momentum transport. In this section, we provide 
a simple prescription for the turbulent stresses T R u ^ b and 
T|T that accounts for the 2D structure of the large-scale 
flow observed in the simulations. 

First, the right hand panel of figure \E\ shows that the Rcf> 
component of the turbulent and viscous stress tensors are 
comparable in amplitude. In addition, both stresses are 
larger by about an order of magnitude than the Z<f> com- 
ponents of the stress Trgjj c and T|^ rfc (compare the vertical 
scales on the two panels of that figure) . This suggests that 
the Z(f> component of the stress is not a fundamental piece 
in the process of setting the flow topology. For example, 
in the derivation of the radial velocity in section 12.31 the 
Rcf> contribution of the viscous stress to the meridional 
variations of the radial velocity writes as 



Ho 
Ro 



R 
Ro 



9+1/2 



3p + 3q + 6 



3q + 9 



(20) 
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Fig. 11. Same as figure [TU1 but for the cases p = —3/2 
(top panel) and p = —5/2 (bottom panel). 



This last expression is very similar to the complete ex- 
pression for vr that is given by Eq. ([T4]) and still indicates 
outward motion in the disk midplane and inward motion 
in the disk upper layers for a wide range of p and q val- 
ues. Obviously, the essence of the meridional circulation 
lies in the R(f> component of the viscous stress. We thus 
adopt, for simplicity, the following prescription for the Z<fi 
component of the turbulent stress: 



Titurb r» 



(21) 



The conservation of angular momentum, Eq. (|15|) . can 
thus be written as 



RpvR 



dR 2c l K d 



dR 



dR 



( r>2rr>turb\ 

>, R 1 r<p ) 



(22) 



where the vertical velocity has been dropped as it is an 
order H/R lower than the radial velocity. Similarly, the 
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vertical variation in f2 appearing on the left hand side of 
the equation above has been neglected because it would 
result in terms of order (H/R) 2 , i.e. much smaller than the 
vertical variations arising as a result of the vertical den- 
sity stratification. Using the variations for the gas density 
given by Eq. ©, Eq. ([2"2"f can then be written as 



vr 
Co 



Ho 
Ra 



R 
Ro 



-p-l/2 p. 



1 



RqPo% 



cxp 



z 2 

2H 2 



(23) 



In the bulk on the disk (\Z\ < 2.5H), the results of the 
simulations (figures [5] and [6]) demonstrate that T R u £ h is 
largely independent of Z. Thus the only dependence with 
Z in the expression above is in the exponential, which ex- 
plains why the sign of the radial velocity does not depend 
on the distance to the midplanc. 

Whether the radial velocity is positive or negative in- 
stead depends on the radial derivative of the turbulent 
stress tensor. If the radial variations of T R u ^ b are steeper 
than R~ 2 , vr is negative everywhere. In the opposite situ- 
ation, it is positive, as found in the simulations presented 
in this paper. For illustrative purposes, we adopt in the 
following a simple power-law radial dependence for T R u £ b 
and writes: 



rpturb 

L R<i> — 



- am c 2 (j£) if \Z\ <2.5H 
otherwise 



(24) 



Su ch a prescription bea rs similarities to the one adopted 
bv lKretke fc Linl (|2010h in their disk models. In Eq. ([2"4")l . 
at is a normalizing factor of the stress that shall not be 
thought of as the standard a parameter. However, on di- 
mensional grounds, at and a are expected to have the 
same order of magnitude. The power-law exponent S is an 
unknown number that depends a priori on the disk prop- 
erties (for example the density and sound speed exponents 
p and q). It is important, at this stage, to stress that such 
a form in the radial variations of the stress tensor has ab- 
solutely no theoretical grounds and should only be viewed 
as a means to pushing the analysis further. Using the ex- 
pression above for the stress tensor, the radial velocity of 
the gas can be written as 



Vr 
Co 



Ho 
Ro 



R 



a-p+i/2 



l^J ex P (^2) (25) 



when \Z\ < 2.5H, while vr vanishes in the disk's upper 
layers. The results of the simulations presented in this 
paper show that vr is positive in the disk midplane, sug- 
gesting S < —2 for the set of disk parameters probed in 
the present paper. It would be tempting at this stage to 
take the actual measured radial variations in the stress 
instead of using the simple power law given by the equa- 
tion above. We chose not to do it for two main reasons: 
first, as detailed in the following sections, these variations 
may be polluted by numerical artifacts and thus may not 
be representative of an actual protoplanetary disk, and 



second, whether a meridional circulation develops does 
not depend on those variations, as shown by Eq. ([23]) . 
Rather, we simply observe that the results of the simula- 
tions indicate a midplane radial velocity on the order of 
5 x 10~ 4 co. In other words, it is around a(Ho/ Rq). This 
suggests in turn that 2{6 + 2)(i?/i? )' 5 ~ p+1/2 is close to 
unity. Eq. (|2"5j) for the radial velocity is compared with the 
simulation results, for the case p = —2, in the left hand 
panel of figure [TUJ In computing vr using that equation, 
we have used a t ~ a ~ 5 x 10~ 3 , H n /R n = 0.1 while the 
remaining terms (with the exception of the exponential!) 
are assumed to be equal to one. Given the simplicity of the 
prescriptions detailed above, the agreement between the 
simulation results and the prediction of Eq. (|25l) is very 
good. 

Finally, the prediction of Eq. (|25[) is further compared 
with the simulation results for the case p = —3/2 (top 
panel) and p = —5/2 (bottom panel) in figure [TO Again, 
the agreement between the simulations and the prediction 
of Eg . ([25]) is satisfactory for the case p = — 5/2, while the 
case p = —3/2 shows larger, but still modest deviations 
from the simple analytical prediction. 

5.2. Comparison with the predictions of standard a 
disk models 

Standard ID a disk theory considers the height and az- 
imuthally integrated version of Eq. ([TUf to describe an- 
gular momentum transport. (All Z-derivatives disappear 
in that procedure.) Using the no tation introduced i n sec- 
tion 13.2.11 it is expressed as ( Papaloizou fc Lir] 1995 : 
Balbus fc Papaloizou|[l999h 



d 

dR 



R 2 Y,T^ C 



IP 



(26) 



Assuming it has the form of a standard viscous stress ten- 
sor, the height and azimuthally integrated R(f> component 
of the stress arc then 



'VISC 

Rcjy 



2tt dR v ' 



Standard a disk models then make the ansatz that it is in 
fact the vertically and azimuthally average kinematic vis- 
cosity V that takes the form given by Eq. (|9]). This is dif- 
ferent from the theory developed in section [2~31 where this 
relation is assumed to hold locally. Using that prescription 
for V in Eq. (|27[) and given the radial dependence of the 
Keplerian angular velocity, the vertically and azimuthally 
averaged stress tensor writes as 



2vr 



T, 



(28) 



For locally isothermal disks such as those considered in the 
present paper, this relation can be equivalently expressed 
in terms of the thermal pressure as 



Ra 



—a / Pd(j>dz . 



(29) 
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This equation shows that the relation between stress and 
pressure holds for vertically integrated quantities in stan- 
dard a disk theory, while the model developed in section [2] 
explicitly assumes that it holds locally: T'^ c — —3/2aP. 
The results presented in this paper have shown that such 
a local relation does not hold between the turbulent stress 
and the thermal pressure, as it would otherwise result 
in meridional circulation developing in turbulent disks. 
Nevertheless, this does not mean that this relation would 
not hold between the vertically averaged turbulent stress 
and the vertically averaged thermal pressure, as is the case 
in standard a disk theory as shown by Eq. (f2l))) . 

A first attempt at answering this question can be made 
by comparing the vertically and azimuthally averaged ra- 
dial velocity predicted, on one side, by ID a disk model 
and, on the other hand, by the numerical simulations 
of protoplanetary disks presented in this paper. In stan- 
dard ID a disk models, the radial profile of the vertically 
integrated radial velocity can be computed by plugging 
Eq. ([29]) into Eq. ® to get 



vr 
Co 



-3a 



Ho 

Rn 
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Ro 



9+1/2 



3q 7 



(30) 



Equation ([3T)|) can also be recovered from Eq. j[T4"]) by a 
simple vertical integration ( Takeuchi fc Lin 20021 ) . For the 
value q = — I considered in the present paper, it vanishes 
and changes sign for p = —2. When p < — 2, Eq. ([50]) pre- 
dicts a positive vertically integrated radial velocity, while 
vr is negative for p > — 2. This behavior is different from 
what is found in the simulations, in which outward bulk 
motions were observed regardless of the value of p. This 
suggests that turbulent disk behave differently than stan- 
dard ID a disk models. 

The systematic outward velocities obtained in the sim- 
ulations can be traced back to a not being a constant 
function of R, but rather systematically decreasing out- 
ward. Indeed, when using the prescription for the turbu- 
lent stress tensor given by Eq. (|2"4"]) . a = TjgT b j p/c^ as 
given by Eq. fTg]) . is found to scale like R s+2 . Constant 
a values thus correspond to disk regions where 5 = — 2. 
In constrast, S < — 2 in those parts of the disks where 
a decreases with R. Using Eq. (|2"5j) . we can thus predict 
outward midplane radial velocities in the region where a 
decreases radially (since S + 2 < there). The tendenct of 
a to decrease with R in our simulations thus explains the 
trend toward finding positive midplane radial velocities. 
As shown by figure [1] in model p = —2, a is almost con- 
stant in the region R S [4.5, 7]. In this region, we thus ex- 
pect to have S ~ 2, and Eq. ([25]) predicts a vanishing radial 
velocity in the disk midplane. This is again consistent with 
the simulation results, as shown in the right hand panel of 
figure 0J Finally, Eq. (|2"5j) predicts that negative midplane 
radial velocities are obtained when S > —2, or equivalently 
when a increases outward. We note that this predictio n is 



indeed reports radially increasing a and inward radial ve- 
locities. 

The systematic positive radial velocities we report here 
might appear suspicious for the long term evolution of 
an accretion disk, so it deserves a comment. We start by 
stressing that a disk models such as those used here also 
exhibit that property for some combination of the expo- 
nents p and q. This is the case for example of the model 
q = — I and p = —2.5, for which Eq. (|3ll|) predicts a 
positive velocity implying outward mass transport. This 
can only be a transient situation. The mass accretion rate 
M = —2-kRY,vr is negative in this case and scales like 
.RP+3(j/2+3_ Thus its exponent is negative, showing that 
the mass accretion rate and its derivative both decrease 
with R in amplitude. On viscous timescales, the radial pro- 
file of S will flatten and accretion will proceed. Ultimately, 
mass will fall onto the star as angular momentum is trans- 
ported outward. Although it is currently impossible to run 
3D MHD simulations on viscous disk timescales, the same 
is most likely true of the turbulent disk that is simulated 
in the present paper. We would thus expect to recover ac- 
cretion if the simulations were run long enough. However, 
this does not prevent a comparison between viscous mod- 
els and simulations when disks are not in steady state, 
because the mean flow in disks is established on a shorter 
timescale than the viscous timescale. The point of the 
present section is precisely to highlight that the predic- 
tions of standard a disk theory differ from the findings of 
numerical simulations. That being said, we caution that 
these simulations only probe such differences for the lim- 
ited time and parameters enabled by current computa- 
tional resources. Whether they are general is beyond the 
scope of this paper and should be the focus of future work. 

We finally comment that the differences we report 
between turbulent and standard a disk model are con- 
sistent with previ o us re sults in the literature. For ex- 
ample, ISano et al.l ([20041 ) report a scaling of the turbu- 
lent stress tensor as P 1 / 4 , and not as P as in standard 
a disk theory . This is consistent with the simulations 
of Lyra et al. ( 20081 ). Along the same lines, Hirose et al 



iave also recently reported the lack of a direct re- 



lation between stress and thermal pressure using shear- 
ing box simulations of turbulent disks. At the same time, 
the results presented in the present paper regarding this 
question should be taken with care. Recent studies have 
shown that the saturation level of the MRI in the un- 
stratified shearing box is strongly affected by explicit dis- 
sipation ( Lesur fc Longarettil 120071 : iFromang et al.l 12007 : 
Simon &: Hawlevll2009 ). In the simulations presented here, 



consistent with the recent results of lFlock et al.l ([201 II ). In 



simulations similar to those presented here, these authors 



all the dissipation is of numerical origin. Thus numerical 
dissipation alone could mask any scaling of the stress with 
the disk parameters. To make things even worse, numeri- 
cal dissipation in these simulations is in principle a func- 
tion of position since the effective radial resolution per disk 
scaleheight varies with radius. Any definite conclusions re- 
garding the precise scaling of the vertically averaged stress 
with pressure should thus wait for better control of dissi- 
pation in global simulations, a difficult task by definition. 
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The present results highlight once more the need for an 
accurate determination of the saturation of the turbulent 
stress as a function of the disk parameters. Such a task is 
beyond the scope of the present paper but should also be 
the focus of future work. 



5.3. Conclusions 

Large-scale meridional circulation has been predicted in 
protoplanetary disks based on 2D viscous disk theory. 
However, the flow in a protoplanetary disk is known to 
be turbulent, most likely because of the MRI. The ques- 
tion raised in this paper is thus simple: does meridional 
circulation exist in turbulent protoplanetary disks? This 
problem has been addressed using a set of numerical simu- 
lations of turbulent protoplanetary disks. The large-scale 
radial flow of the gas was computed by averaging the re- 
sults in the azimuthal direction, as well as on long time 
intervals. The results are found to disagree with 2D viscous 
disk theory. There is no sign of any meridional circulation 
in the disk. Instead, the sign of the radial velocity is found 
to be constant in the bulk of the disk (\Z\ < 2.5H). This 
is shown to come from the vertical behavior of the turbu- 
lent stress tensor. Instead of scaling like the local thermal 
pressure as assumed in 2D viscous disk theory, the tur- 
bulent stress is constant around the disk midplane before 
dropping in the disk corona, where thermal and magnetic 
pressure become comparable. Such a vertical profile of the 
turbulent stress has been found already in many numeri- 
cal simulations, regarless of their local and global nature. 
It is thus a robust property of the flow in turbulent disks. 
In light of these properties of the stress, the conclusion 
that no meridional circulation exists in fully ionized and 
turbulent protoplanetary disks appears unavoidable. 

This has important consequences for the trans- 
port of solids in the solar nebula. The relevance of 
models in which crystalline dust radial transport is 
largely based on meridiona l circulation can be ques- 
tioned iKeiler fc Gail! l2004t iTscharnuter fc Gail 12007 : 



Cieslall2007l.l2009l iHughes fc ArmitagefeOlOh . In addition 
the simulations also raise the possibility of an outward ra- 
dial mass flux unexpected from a standard a disk model. 
This opens up the possibility of an efficient mechanism 
to transport solid particles outward. Given the current 
state of the art of global numerical simulations of pro- 
toplanetary disks, the viability of such a feature of the 
flow nonetheless remains very uncertain. Future studies of 
MRI-powered MHD turbulence should aim at better con- 
straining the nature of the large scale flow in disks and 
the saturation amplitude of the MRI as a function of the 
disk's properties (surface densities, temperature). 

It is also important to realize that the simulations pre- 
sented here do not rule out the possibility of a large- 
scale outward flow in protoplanetary disks. They only 
demonstrate it does not exist in fully turbulent disks 
and that a viscous modeling of such disks leads to er- 
roneous results. There are a number of possibilities by 



which a meridional flow might develop in real proto- 
planetary disks. First, protoplanetary disks are known to 
be largely neutral and most likely harbor a dead zone 
around their equatorial plane in which the flow is sta- 
ble to the MRI. While the exact size and dynamical sta- 



both on chemical (Fromang et alJ 2002: Ilgner & Nelson! 


2006bllal) and dynamical effects (F 


eming & Stone! 20031: 


Turner et al. 2007: Ilgner & Nelson 


20081). its verv exis- 



tence in the planet -forming region seems difficult to avoid. 
The large-scale flow in such disks has never been con- 
sidered and might lead to unexpected results that could 
affect t he large-scal e trans port of solid particles. For ex- 
ample, Bai fc Stone! ( 2010h recently showed that the fast 



inward drift of large particles causes an outflow of gas 
in the dead zone midplane that could carry small solids 
like CAIs with it. Finally, one should bear in mind the 
idealized nature and the limited sampling of the param- 
eter space of the simulations presented here. Additional 
physical ingredients (ambipolar diffusion, hall effects) or 
different disk properties (temperature radial and vertical 
profile, magnetic field configuration) might all change the 
vertical profile of the stress tensor and create a large scale 
outflow in the bulk of the disk. 
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Appendix A: Off-centered Keplerian viscous ring 
spreading in spherical geometry 

o To test the implementation of the viscous stress ten- 
sor in our purely hydrodynamical codes, we have designed 
a set up, in spherical coordinates (r,0,<j>), in which the 
two components of the stress tensor that potentially con- 
tribute to meridional circulation are within the same order 
of magnitude. It consists of a disk orbiting within a fixed 
potential $ that has the following form: 



GM 



rsmf 



A(r cos ( 



(A.l) 



where Zq and A are constants. This potential is Keplerian 
along the cylindrical radius R = r sin 9, and harmonic in 
the vertical direction Z = r cos 9, with a potential min- 
imum at Zq 7^ 0. An equilibrium setup for this poten- 
tial can therefore correspond, for a globally isothermal 
gas with sound speed c s , to an approximately Keplerian 
disk with its equator lying at Z = Z$ and with a uniform 
thickness given by H — c s /u> z , where uj 2 = d 2 Q/dz 2 = 2A. 
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Namely, we initialize our disk as 



p{r,6,<t>) 



-3/4 



(A.2) 



exp 



- l) 2 (Z - Z 



TO 



2H 2 



where x = R/Rq and tq = \2vtjR\, which is the ra- 
tio of the time and of the v iscous timescale of the disk 
( Lvnden-Bell fe Pringld 1974 ) . t§ is set to be much lower 
than unity in the expression above. Equation (|A.2[) cor- 
responds to a narrow Gaussian ring of mass to, centered 
on i?o, in vertical hydrostatic equilibrium, so that it also 
has a Gaussian profile, centered on Zq, in the vertical 
di rection. Eq. (|A.2I) stems from the classical expression 
of iLvnden-Bell fc Pringld (|l974l) . which describes the ra- 
dial spread of an initially infinitely narrow viscous ring, 
in which we approximate the Bessel function J 1 / 4 (2a;/To) 
with exp(2a;/ro)/ ^4ttx/tq, using the fact that its argu- 
ment 2x/to is large. This means that we start with a ring 
that has already undergone some viscous radial spread, al- 
though for a time short compared to its viscous timescale, 
so that its width is small compared to its radius. A 
meridional cut of the initial density field is represented 
in Figure EU 

We initialize the azimuthal velocity as follows: 

1/2 



r gm. 



2(x - l)x 



To 



(A.3) 



where the last two terms correspond to the rotational sup- 
port provided by the radial pressure gradient, and we ini- 
tialize the cylindrical radial velocity as 

'3 6 



VR 



R 



H x(x - 1) 



TO 



(A.4) 



from which we set v r = VRSm.9 and vg = vrcosO. 
Equation (|A.4[) stems from Eq. (fT3"f , in which the Z deriva- 
tive is set to 0, as no vertical dependence of il is expected 
in this setup, each layer at a given altitude being a replica 
of the equatorial layer, with a density ratio that depends 
on the altitude. 

For our fiducial calculation, we adopt c s = 
lO-'^^GM/fio) 1 / 2 , A = 10- 2 /3(GM/i?3), hence H = 
V0.15i?o, and we start with a narrow ring correspond- 
ing to To = 0.018, with an equator at Zq = 1. Our mesh 
contains 50 x 50 zones, regularly spaced in radius from 
R = 0.8i?o to R — 2R , and regularly spaced in colatitude 
from 9 = 0.35 to 9 = 1.3. Since the azimuthal velocities 
are taken into account, our setup is "2D1/2" in nature. 

We present in Figure IA.2I the results of this calcu- 
lation, which show that the radial spread of the ring is 
correctly captured on the spherical mesh. We performed 
subsidiary calculations in which we manually cancel out 
one of the two stress components that have an impact on 
the torque felt by a ring of material (hence on the rate 
at which the material spreads), i.e. either Tg^ c or T™ c . 
That the resulting curves of these two additional calcula- 
tions display a significantly different result than the cor- 
rect result and that they yield a similar peak value at 



1.5 



l.O 



0.5 




0.5 



1.0 



1.5 



2.0 



Fig. A.l. Initial density field (gray levels) superimposed 
on the spherical mesh used in the calculation. The dashed 
line represents the "equator" (i.e. the potential minimum 
in the Z direction), while the thick arrows depict the ex- 
pected trend of the ring's material to spread radially. 



the same date shows that these two components of the 
stress tensor play roles of comparable magnitude on the 
spread of the ring. This can be understood as the spheri- 
cal coordinate system, in the vicinity of the center of the 
ring at Z = Zq, is approximately tilted by 45° with re- 
spect to the equator. We finally note th at the original 
derivation of ILvnden-Bell fc Pringle (Il974h does not con- 
template pressure effects. Those are relatively substantial 
in the disk considered here, and might account for the 
tiny discrepancy between the simulated and expected pro- 
files in Figure lA~2l except near the boundaries, where the 
discrepancy likely arises from the accumulation of mate- 
rial due to the reflexive boundary conditions that we use, 
which prevent material from flowing out of the grid. 
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